home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / AA_51.ZIP / ONEPTUNE.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  3KB  |  122 lines
  1. /* Orbital elements and perturbations for the planet Neptune
  2.  * using formulas given by Meeus
  3.  */
  4. #include "planet.h"
  5. #include "kep.h"
  6.  
  7. static double sinG = 0.0;
  8. static double cosG = 0.0;
  9.  
  10. int pjup();
  11.  
  12. int oneptune(e,J)
  13. struct orbit *e;
  14. double J;
  15. {
  16. double p, q;
  17.  
  18. e->epoch = J;
  19. manoms(J);
  20. pjup();
  21.  
  22. f = (-0.000070*T + 218.46134)*T + 37.73063;
  23. M8 = mod360(f);
  24. e->M = M8;
  25. e->a = 30.10957;
  26. e->ecc = ( -0.000000002*T + 0.000006330)*T + 0.00899704;
  27. #if OFDATE
  28. e->equinox = J;
  29. e->i = ( -0.0000091*T - 0.0095436)*T + 1.779242;
  30. f = (( 0.000004113*T + 0.00014095)*T + 0.3256394)*T + 276.045975;
  31. e->w = mod360(f);
  32. f = (( -0.000004718*T + 0.00024987)*T + 1.0989350)*T + 130.681389;
  33. e->W = mod360(f);
  34. f = (( -0.00000060*T + 0.0003205)*T + 219.885914)*T + 84.457994;
  35. e->L = mod360(f);
  36. #else
  37. e->equinox = J2000;
  38. e->i = ((1.8e-8*T - 2.27e-6)*T - 0.0000144)*T + 1.769715;
  39. f = ((2.226e-6*T + 3.849e-5)*T + 0.0368127)*T + 276.335328;
  40. e->w = mod360(f);
  41. f = (( -2.858e-6*T + 4.428e-5)*T - 0.0084187)*T + 131.788486;
  42. e->W = mod360(f);
  43. f = e->M + e->w + e->W;
  44. e->L = mod360(f);
  45. #endif
  46.  
  47. G = 83.76922 + 218.4901*T;
  48. G = mod360(G)*DTR;
  49. H = modtp(2.0*G - S);
  50. ze = modtp(G - P);
  51. eta = modtp(G - Q);
  52. th = modtp(G - S);
  53. sinth = sin(th);
  54. costh = cos(th);
  55. sin2th = 2.0*sinth*costh;
  56. cos2th = costh*costh - sinth*sinth;
  57. sin1H = sin(H);
  58. cos1H = cos(H);
  59. sin2H = 2.0*sin1H*cos1H;
  60. cos2H = cos1H*cos1H - sin1H*sin1H;
  61. sinG = sin(G);
  62. cosG = cos(G);
  63.  
  64. /* perturbations in mean longitude */        
  65. p = (-0.589833 + 0.001089*nu)*sin1H;
  66. p = p +(-0.056094 + 0.004658*nu)*cos1H;
  67. p = p -0.024286*sin2H;
  68. e->L += p;
  69.  
  70. /* perturbations in the perihelion */
  71. q = 0.024039*sin1H
  72.    +0.006206*sin2H
  73.    -0.025303*cos1H
  74.    -0.005992*cos2H;
  75. e->w += q;
  76. e->M = M8 + p - q/(e->ecc);       
  77.  
  78. /* perturbations in eccentricity */
  79. q = .0004389*sin1H
  80.    +.0001129*sin2H
  81.    +.0004262*cos1H
  82.    +.0001089*cos2H;
  83. e->ecc += q;
  84.  
  85. /* correction to the semimajor axis */
  86. q = -0.000817*sin1H
  87.     +0.008189*cos1H
  88.     +0.000781*cos2H;
  89. e->a += q;        
  90. return(0);
  91. }
  92.  
  93.  
  94.  
  95. int cneptune(e)
  96. struct orbit *e;
  97. {
  98. double q;
  99.  
  100. /* correction to the true longitude */
  101. q = -0.009556*sin(ze)
  102.     -0.005178*sin(eta)
  103.     +0.002572*sin2th;
  104. q = q -0.002972*cos2th*sinG;
  105. q = q -0.002833*sin2th*cosG;
  106. e->L += q*DTR;
  107.  
  108. /* correction to the heliocentric latitude */
  109. q = 0.000336*cos2th*sinG;
  110. q = q + 0.000364*sin2th*cosG;
  111. e->plat = q * DTR;
  112.  
  113. /* correction to the radius vector */
  114. q = -0.040596
  115.     +0.004992*cos(ze)
  116.     +0.002744*cos(eta)
  117.     +0.002044*costh
  118.     +0.001051*cos2th;
  119. e->r += q;
  120. return(0);
  121. }
  122.